library(INLA)
library(rgdal)
library(ggplot2)
library(viridis)
library(gridExtra)
library(raster)
library(ggspatial)
library(rgeos)
# Xylella fastidiosa occurrence data in the demarcated area of Alicante (Spain)
data_Xf <- read.table("./Data/data_Xf.txt")
# Coordinate reference system
newproj <- "+proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs"
## -- Boundaries -- ##
# Study area
bound <- readOGR('./Data/Boundaries', 'bound')
# Mountain barrier boundary
MB <- readOGR('./Data/Boundaries', 'bound_MB')
bound_MB <- gDifference(bound, MB, checkValidity = 2L)
# Continuous barrier boundary
CB <- readOGR('./Data/Boundaries', 'bound_CB')
bound_CB <- gDifference(bound, CB, checkValidity = 2L)
# Discontinuous barrier boundary
DB <- readOGR('./Data/Boundaries', 'bound_DB')
bound_DB <- gDifference(bound, DB, checkValidity = 2L)
## -- Remove data inside the barrier -- ##
# Mountain barrier
out.barrier = over(bound_MB,
SpatialPoints(cbind(data_Xf$X, data_Xf$Y), proj4string = CRS(newproj)),
returnList = T)[[1]]
data_M <- data_Xf[out.barrier,]
# write.table(data_M, './Data/data_Xf_M.txt')
# Perimeter barriers
out.barrier = over(bound_CB,
SpatialPoints(cbind(data_Xf$X, data_Xf$Y), proj4string = CRS(newproj)),
returnList = T)[[1]]
data_B <- data_Xf[out.barrier,]
# write.table(data_B, './Data/data_Xf_B.txt')
## -- Read data -- ##
# Mountain barrier data
data_M <- read.table("./Data/data_Xf_M.txt")
# Perimeter barriers data
data_B <- read.table("./Data/data_Xf_B.txt")
## -- Mesh -- ##
# Stationary
mesh_ST <- inla.mesh.2d(
boundary = bound,
max.edge = c(750, 3500),
offset = c(0.1,-0.1),
cutoff = 750
)
# Mountain barrier
mesh_MB <- inla.mesh.2d(
boundary = bound_MB,
max.edge = c(750, 3500),
offset = c(0.1, -0.1),
cutoff = 750
)
# Continuous barrier
mesh_CB <- inla.mesh.2d(
boundary = bound_CB,
max.edge = c(750, 3500),
offset = c(0.1, -0.1),
cutoff = 750
)
# Discontinuous barrier
mesh_DB <- inla.mesh.2d(
boundary = bound_DB,
max.edge = c(750, 3500),
offset = c(0.1, -0.1),
cutoff = 750
)
data <- data_Xf
mesh <- mesh_ST
# SPDE model definition
size <- min(c(diff(range(data$X)), diff(range(data$Y))))
range0 <- size / 2
spde <- inla.spde2.pcmatern(
# Mesh
mesh = mesh,
# P(practic.range < range0) = 0.5
prior.range = c(range0, 0.5),
# P(sigma > 10) = 0.01
prior.sigma = c(10, 0.01)
)
# Projector matrix
A.est <- inla.spde.make.A(mesh, loc = cbind(data$X, data$Y))
# Data stack
stk.est <- inla.stack(
data = list(y = data$ResLab),
A = list(A.est, 1),
effects = list(spatial = 1:mesh$n,
beta0 = rep(1, nrow(data))),
tag = 'est'
)
# Model fitting
formula.1 <- y ~ -1 + beta0 + f(spatial, model = spde)
model.est_ST <- inla(
formula.1,
data = inla.stack.data(stk.est),
family = "binomial" ,
control.compute = list(
dic = TRUE,
cpo = TRUE,
waic = TRUE,
return.marginals = TRUE
),
control.predictor = list(A = inla.stack.A(stk.est),
compute = TRUE),
num.threads = 2,
verbose = FALSE
)
saveRDS(model.est_ST, "./Models/ST_model_est.rds")
# Read saved model
model.est_ST <- readRDS("./Models/ST_model_est.rds")
# Projection on a grid
dxy <- apply(bbox(bound), 1, diff)
r <- dxy[1] / dxy[2]
m <- 120
proj.grid.mat <-
inla.mesh.projector(
mesh,
xlim = bbox(bound)[1, ],
ylim = bbox(bound)[2, ],
dims = c(r, 1) * m
)
# NA to the values outside boundary
ov <-
over(
SpatialPoints(proj.grid.mat$lattice$loc, bound@proj4string),
SpatialPolygons(bound@polygons, proj4string = bound@proj4string)
)
i.temp <- is.na(ov)
# Projector matrix to predict
A.pred <-
inla.spde.make.A(mesh, loc = proj.grid.mat$lattice$loc[!i.temp, ])
# Stack to predict
stk.pred <- inla.stack(
data = list(y = NA),
A = list(A.pred, 1),
effects = list(spatial = 1:mesh$n,
beta0 = rep(1, dim(A.pred)[1])),
tag = 'pred'
)
stk <- inla.stack(stk.est, stk.pred)
# Prediction
model.pred_ST <- inla(
formula.1,
data = inla.stack.data(stk),
family = "binomial",
control.predictor = list(
A = inla.stack.A(stk),
compute = TRUE,
link = 1
),
control.mode = list(theta = model.est_ST$mode$theta,
restart = TRUE),
num.threads = 3,
verbose = TRUE
)
saveRDS(model.pred_ST, "./Models/ST_model_pred.rds")
# Read saved model
model.pred_ST <- readRDS("./Models/ST_model_pred.rds")
# Index of the predicted values
idx <- inla.stack.index(stk, 'pred')$data
# Matrix to visualize
prob.mean_ST <-
prob.sd_ST <- matrix(NA, proj.grid.mat$lattice$dims[1],
proj.grid.mat$lattice$dims[2])
prob.mean_ST[!i.temp] <-
c(model.pred_ST$summary.fitted.values$mean[idx])
prob.sd_ST[!i.temp] <-
c(model.pred_ST$summary.fitted.values$sd[idx])
summary(model.est_ST)$fixed
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## beta0 -1.68 0.245 -2.205 -1.666 -1.233 -1.641 0
Â
summary(model.est_ST)$hyperpar
Â
Â
# Raster of the posterior predictive mean for the difference between the models
mean_pred_ST <- matrix_to_raster(prob.mean_ST, proj.grid.mat = proj.grid.mat)
mean_pred_ST<-mask(mean_pred_ST, bound)
data <- data_M
bound_barriers <- bound_MB
mesh <- mesh_MB
# Polygon with the barrier
in.tri <- inla.over_sp_mesh(bound_barriers,
y = mesh,
type = "centroid",
ignore.CRS = TRUE)
num.tri <- length(mesh$graph$tv[, 1])
barrier.tri <- setdiff(1:num.tri, in.tri)
poly.barrier <- inla.barrier.polygon(mesh,
barrier.triangles = barrier.tri)
# SPDE model definition
size <- min(c(diff(range(data$X)), diff(range(data$Y))))
range0 <- size / 2
barrier.model <- inla.barrier.pcmatern(
mesh,
barrier.triangles =
barrier.tri,
prior.range = c(range0, 0.5),
prior.sigma = c(10, 0.01)
)
# Projector matrix
A.est <- inla.spde.make.A(mesh, loc = cbind(data$X, data$Y))
# Data stack
stk.est <- inla.stack(
data = list(y = data$ResLab),
A = list(A.est, 1),
effects = list(spatial = 1:mesh$n,
beta0 = rep(1, nrow(data))),
tag = 'est'
)
# Model fitting
formula.1 <- y ~ -1 + beta0 + f(spatial, model = barrier.model)
model.est_MB <- inla(
formula.1,
data = inla.stack.data(stk.est),
family = "binomial" ,
control.compute = list(
dic = TRUE,
cpo = TRUE,
waic = TRUE,
return.marginals = TRUE
),
control.predictor = list(A = inla.stack.A(stk.est),
compute = TRUE),
num.threads = 2,
verbose = FALSE
)
saveRDS(model.est_MB, "./Models/MB_model_est.rds")
# Read saved model
model.est_MB <- readRDS("./Models/MB_model_est.rds")
# Projection on a grid
dxy <- apply(bbox(bound_barriers), 1, diff)
r <- dxy[1] / dxy[2]
m <- 120
proj.grid.mat <-
inla.mesh.projector(
mesh,
xlim = bbox(bound_barriers)[1, ],
ylim = bbox(bound_barriers)[2, ],
dims = c(r, 1) * m
)
# NA to the values outside boundary
ov <-
over(
SpatialPoints(proj.grid.mat$lattice$loc, bound@proj4string),
SpatialPolygons(bound@polygons, proj4string = bound@proj4string)
)
i.temp <- is.na(ov)
# Projector matrix to predict
A.pred <-
inla.spde.make.A(mesh, loc = proj.grid.mat$lattice$loc[!i.temp, ])
# Stack to predict
stk.pred <- inla.stack(
data = list(y = NA),
A = list(A.pred, 1),
effects = list(spatial = 1:mesh$n,
beta0 = rep(1, dim(A.pred)[1])),
tag = 'pred'
)
stk <- inla.stack(stk.est, stk.pred)
# Prediction
model.pred_MB <- inla(
formula.1,
data = inla.stack.data(stk),
family = "binomial",
control.predictor = list(
A = inla.stack.A(stk),
compute = TRUE,
link = 1
),
control.mode = list(theta = model.est_MB$mode$theta,
restart = TRUE),
num.threads = 3,
verbose = TRUE
)
saveRDS(model.pred_MB, "./Models/MB_model_pred.rds")
# Read saved model
model.pred_MB <- readRDS("./Models/MB_model_pred.rds")
# Index of the predicted values
idx <- inla.stack.index(stk, 'pred')$data
# Matrix to visualize
prob.mean_MB <-
prob.sd_MB <- matrix(NA, proj.grid.mat$lattice$dims[1],
proj.grid.mat$lattice$dims[2])
prob.mean_MB[!i.temp] <-
c(model.pred_MB$summary.fitted.values$mean[idx])
prob.sd_MB[!i.temp] <-
c(model.pred_MB$summary.fitted.values$sd[idx])
summary(model.est_MB)$fixed
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## beta0 -1.612 0.227 -2.093 -1.601 -1.193 -1.581 0
Â
summary(model.est_MB)$hyperpar
Â
Standard deviation of the spatial effect \(= e^{\theta_1}\)
Range \(= e^{\theta_2}\)
pander::pander(exp(model.est_MB$summary.hyperpar[,c(1,3,4,5)]), row.names=c("SD", "Range"))
| Â | mean | 0.025quant | 0.5quant | 0.975quant |
|---|---|---|---|---|
| SD | 1.433 | 1.205 | 1.432 | 1.71 |
| Range | 3861 | 2919 | 3844 | 5212 |
Â
Â
# Raster of the posterior predictive mean for the difference between the models.
mean_pred_MB <-
matrix_to_raster(prob.mean_MB, proj.grid.mat = proj.grid.mat)
mean_pred_MB <- mask(mean_pred_MB, bound)
data <- data_B
bound_barriers <- bound_CB
mesh <- mesh_CB
# Polygon with the barrier
in.tri <- inla.over_sp_mesh(bound_barriers,
y = mesh,
type = "centroid",
ignore.CRS = TRUE)
num.tri <- length(mesh$graph$tv[, 1])
barrier.tri <- setdiff(1:num.tri, in.tri)
poly.barrier <- inla.barrier.polygon(mesh,
barrier.triangles = barrier.tri)
# SPDE model definition
size <- min(c(diff(range(data$X)), diff(range(data$Y))))
range0 <- size / 2
barrier.model <- inla.barrier.pcmatern(
mesh,
barrier.triangles =
barrier.tri,
prior.range = c(range0, 0.5),
prior.sigma = c(10, 0.01)
)
# Projector matrix
A.est <- inla.spde.make.A(mesh, loc = cbind(data$X, data$Y))
# Data stack
stk.est <- inla.stack(
data = list(y = data$ResLab),
A = list(A.est, 1),
effects = list(spatial = 1:mesh$n,
beta0 = rep(1, nrow(data))),
tag = 'est'
)
# Model fitting
formula.1 <- y ~ -1 + beta0 + f(spatial, model = barrier.model)
model.est_CB <- inla(
formula.1,
data = inla.stack.data(stk.est),
family = "binomial" ,
control.compute = list(
dic = TRUE,
cpo = TRUE,
waic = TRUE,
return.marginals = TRUE
),
control.predictor = list(A = inla.stack.A(stk.est),
compute = TRUE),
num.threads = 2,
verbose = FALSE
)
saveRDS(model.est_CB, "./Models/CB_model_est.rds")
# Read saved model
model.est_CB <- readRDS("./Models/CB_model_est.rds")
# Projection on a grid
dxy <- apply(bbox(bound_barriers), 1, diff)
r <- dxy[1] / dxy[2]
m <- 120
proj.grid.mat <-
inla.mesh.projector(
mesh,
xlim = bbox(bound_barriers)[1, ],
ylim = bbox(bound_barriers)[2, ],
dims = c(r, 1) * m
)
# NA to the values outside boundary
ov <-
over(
SpatialPoints(proj.grid.mat$lattice$loc, bound@proj4string),
SpatialPolygons(bound@polygons, proj4string = bound@proj4string)
)
i.temp <- is.na(ov)
# Projector matrix to predict
A.pred <-
inla.spde.make.A(mesh, loc = proj.grid.mat$lattice$loc[!i.temp, ])
# Stack to predict
stk.pred <- inla.stack(
data = list(y = NA),
A = list(A.pred, 1),
effects = list(spatial = 1:mesh$n,
beta0 = rep(1, dim(A.pred)[1])),
tag = 'pred'
)
stk <- inla.stack(stk.est, stk.pred)
# Prediction
model.pred_CB <- inla(
formula.1,
data = inla.stack.data(stk),
family = "binomial",
control.predictor = list(
A = inla.stack.A(stk),
compute = TRUE,
link = 1
),
control.mode = list(theta = model.est_CB$mode$theta,
restart = TRUE),
num.threads = 3,
verbose = TRUE
)
saveRDS(model.pred_CB, "./Models/CB_model_pred.rds")
# Read saved model
model.pred_CB <- readRDS("./Models/CB_model_pred.rds")
# Index of the predicted values
idx <- inla.stack.index(stk, 'pred')$data
# Matrix to visualize
prob.mean_CB <-
prob.sd_CB <- matrix(NA, proj.grid.mat$lattice$dims[1],
proj.grid.mat$lattice$dims[2])
prob.mean_CB[!i.temp] <-
c(model.pred_CB$summary.fitted.values$mean[idx])
prob.sd_CB[!i.temp] <-
c(model.pred_CB$summary.fitted.values$sd[idx])
summary(model.est_CB)$fixed
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## beta0 -1.79 0.375 -2.633 -1.755 -1.149 -1.691 0
Â
summary(model.est_CB)$hyperpar
Â
Standard deviation of the spatial effect \(= e^{\theta_1}\)
Range \(= e^{\theta_2}\)
pander::pander(exp(model.est_CB$summary.hyperpar[,c(1,3,4,5)]), row.names=c("SD", "Range"))
| Â | mean | 0.025quant | 0.5quant | 0.975quant |
|---|---|---|---|---|
| SD | 1.496 | 1.202 | 1.494 | 1.875 |
| Range | 6141 | 4296 | 6103 | 9043 |
Â
# Raster of the posterior predictive mean for the difference between the models.
mean_pred_CB <-
matrix_to_raster(prob.mean_CB, proj.grid.mat = proj.grid.mat)
mean_pred_CB <- mask(mean_pred_CB, bound)
data <- data_B
bound_barriers <- bound_DB
mesh <- mesh_DB
# Polygon with the barrier
in.tri <- inla.over_sp_mesh(bound_barriers,
y = mesh,
type = "centroid",
ignore.CRS = TRUE)
num.tri <- length(mesh$graph$tv[, 1])
barrier.tri <- setdiff(1:num.tri, in.tri)
poly.barrier <- inla.barrier.polygon(mesh,
barrier.triangles = barrier.tri)
# SPDE model definition
size <- min(c(diff(range(data$X)), diff(range(data$Y))))
range0 <- size / 2
barrier.model <- inla.barrier.pcmatern(
mesh,
barrier.triangles =
barrier.tri,
prior.range = c(range0, 0.5),
prior.sigma = c(10, 0.01)
)
# Projector matrix
A.est <- inla.spde.make.A(mesh, loc = cbind(data$X, data$Y))
# Data stack
stk.est <- inla.stack(
data = list(y = data$ResLab),
A = list(A.est, 1),
effects = list(spatial = 1:mesh$n,
beta0 = rep(1, nrow(data))),
tag = 'est'
)
# Model fitting
formula.1 <- y ~ -1 + beta0 + f(spatial, model = barrier.model)
model.est_DB <- inla(
formula.1,
data = inla.stack.data(stk.est),
family = "binomial" ,
control.compute = list(
dic = TRUE,
cpo = TRUE,
waic = TRUE,
return.marginals = TRUE
),
control.predictor = list(A = inla.stack.A(stk.est),
compute = TRUE),
num.threads = 2,
verbose = FALSE
)
saveRDS(model.est_DB, "./Models/DB_model_est.rds")
# Read saved model
model.est_DB <- readRDS("./Models/DB_model_est.rds")
# Projection on a grid
dxy <- apply(bbox(bound_barriers), 1, diff)
r <- dxy[1] / dxy[2]
m <- 120
proj.grid.mat <-
inla.mesh.projector(
mesh,
xlim = bbox(bound_barriers)[1, ],
ylim = bbox(bound_barriers)[2, ],
dims = c(r, 1) * m
)
# NA to the values outside boundary
ov <-
over(
SpatialPoints(proj.grid.mat$lattice$loc, bound@proj4string),
SpatialPolygons(bound@polygons, proj4string = bound@proj4string)
)
i.temp <- is.na(ov)
# Projector matrix to predict
A.pred <-
inla.spde.make.A(mesh, loc = proj.grid.mat$lattice$loc[!i.temp, ])
# Stack to predict
stk.pred <- inla.stack(
data = list(y = NA),
A = list(A.pred, 1),
effects = list(spatial = 1:mesh$n,
beta0 = rep(1, dim(A.pred)[1])),
tag = 'pred'
)
stk <- inla.stack(stk.est, stk.pred)
# Prediction
model.pred_DB <- inla(
formula.1,
data = inla.stack.data(stk),
family = "binomial",
control.predictor = list(
A = inla.stack.A(stk),
compute = TRUE,
link = 1
),
control.mode = list(theta = model.est_DB$mode$theta,
restart = TRUE),
num.threads = 3,
verbose = TRUE
)
saveRDS(model.pred_DB, "./Models/DB_model_pred.rds")
# Read saved model
model.pred_DB <- readRDS("./Models/DB_model_pred.rds")
# Index of the predicted values
idx <- inla.stack.index(stk, 'pred')$data
# Matrix to visualize
prob.mean_DB <-
prob.sd_DB <- matrix(NA, proj.grid.mat$lattice$dims[1],
proj.grid.mat$lattice$dims[2])
prob.mean_DB[!i.temp] <-
c(model.pred_DB$summary.fitted.values$mean[idx])
prob.sd_DB[!i.temp] <-
c(model.pred_DB$summary.fitted.values$sd[idx])
summary(model.est_CB)$fixed
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## beta0 -1.79 0.375 -2.633 -1.755 -1.149 -1.691 0
Â
summary(model.est_CB)$hyperpar
Â
Standard deviation of the spatial effect \(= e^{\theta_1}\)
Range \(= e^{\theta_2}\)
pander::pander(exp(model.est_CB$summary.hyperpar[,c(1,3,4,5)]), row.names=c("SD", "Range"))
| Â | mean | 0.025quant | 0.5quant | 0.975quant |
|---|---|---|---|---|
| SD | 1.496 | 1.202 | 1.494 | 1.875 |
| Range | 6141 | 4296 | 6103 | 9043 |
Â
# Raster of the posterior predictive mean for the difference between the models.
mean_pred_DB <-
matrix_to_raster(prob.mean_DB, proj.grid.mat = proj.grid.mat)
mean_pred_DB <- mask(mean_pred_DB, bound)